%Spin Echo SS-PARSE
if rep == 0
    WW64=(sqrt(XX.^2+XX'.^2)<fov/2);
    reso=64;
    fov=12.8;
    JR=zeros(1,200);de=JR;CC=JR;
    del=4e-09; 
    psis=zeros(NLS,NNN); 
    fr=zeros(64,64);
    amp=fr;
    expol = zeros(64,64);
    expt = zeros(64,64);
    twopi=2*pi;
    for nl=1:NLS
        psis(nl,:)=exp(coex*(kx(1:NNN)*xa(nl)+ky(1:NNN)*ya(nl)));
    end
end 
rep = 1;
increml=[2.e-06;1e-06;6e-07;8.e-08;3e-08];nincr=1;
wfr = 0.4;
wfi = 0.8;
wfr1 = 0.3;
NLIST=[(1:10).'*1188;13000;15000;17000;19000;21000;23763];
FLIST =[0.1;0.06;0.06;0.05;0.05;0.05;0.04;0.04;0.03;0.03;0.02;0.01;0.01;0.005;0.005;0.004];     
NLENGTHS=size(NLIST,1);
NIT=120;
R2f=zeros(64,64,NIT);
R2P = zeros(64,64,NIT);
 p0((1:NLS)+NLS)=real(ptem(1:NLS))/wfr+1i*imag(ptem(1:NLS))/wfi;
p01(1:NLS) = newtemp(1:NLS)/wfr1;
p=p0;
p2 = p01;
for nl=1:NLS
    if abs(wfr*real(p(nl+NLS))/delt) > 150;p(nl)=0.001*p(nl);p(nl+NLS)=0;end
    if real(p(nl+NLS))>0; p(nl+NLS)=0+1i*imag(p(nl+NLS));end
end
for nl = 1:NLS
    if abs(wfr1*p2(nl)/delt)>150;
        p2(nl)=0.001*p2(nl);
    end
    if p2(nl)>0
        p2(nl) = 0;
    end
end
% WR=sqrt(XX.^2+XX'.^2);
WRA=WW64;
%WRX=WW64.*(.54+.46*cos(pi*(WR./6.4)));
% WRA=WRX;

FID3=FIDD;
NNN0=0;niter=1;CCold=1e+15;CCnew=1e+14;
tic
for nlen=1:NLENGTHS
    NNN1=NLIST(nlen);frac=FLIST(nlen);
    while abs(CCold-CCnew) > frac*CCold || NNN1 ~= NNN0
        p1=p;
        p3=p2;
        % DJ evaluation
        DJ=zeros(1,5*NLS);
        sesto=xyvalcs(p1,p3,psis,NNN1,NLS,wfr,wfi,wfr1,time1,time2);
        sigdiff=(FID3(1:NNN1)-sesto);
        sigdif0 = sigdiff;
        J0=sigdif0*sigdif0';
        nincr=1;
        if NNN1 > 2000
            nincr=2;
        end
        if NNN1 >3600
            nincr=3;
        end
        if NNN1 >6000
            nincr=4;
        end
        if NNN1 >8000
            nincr=5;
        end
        increm=increml(nincr);
        at=p(1:NLS).';
        bt=wfr*real(p(NLS+1:2*NLS)).'+1i*wfi*imag(p(NLS+1:2*NLS)).';
        ct = wfr1*p2(1:NLS).';
        clde = exp(wfr*abs(time2(1:NNN1))*del);
        clide = exp(wfi*1i*(time2(1:NNN1))*del);
        clide1 = exp(wfr1*abs(time1(1:NNN1))*del);
        figure(3);subplot(2,1,1);title('computing gradients');drawnow;
        for nl = 1: NLS 
            clz = exp((abs(time2(1:NNN1))*real(bt(nl)))+(1i*imag(bt(nl))*time2(1:NNN1)));
            % clz = exp(bt(nl)*abs(time2(1:NNN1)));
            c10 = clz.*psis(nl,1:NNN1);
            cnew =  exp(abs(time1(1:NNN1))*ct(nl));          

            sigdif1 = sigdiff+at(nl).*c10 .*cnew;
            sesub0=(at(nl)+del)*c10.*cnew;
            sigdif2=(sigdif1-sesub0);
            DJ(nl)=(sigdif2*sigdif2'-J0)/del;

            sesub0=(at(nl)+1i*del)*c10.*cnew;
            sigdif2 = (sigdif1-sesub0);
            DJ(nl+2*NLS)=(sigdif2*sigdif2'-J0)/del;

            sesub0=at(nl)*(c10.*clde).*cnew;
            sigdif2 = (sigdif1-sesub0);
            DJ(nl+NLS)=(sigdif2*sigdif2'-J0)/del;

            sesub0=at(nl)*(c10.*clide).*cnew;
            sigdif2 = (sigdif1-sesub0);
            DJ(nl+3*NLS)=(sigdif2*sigdif2'-J0)/del;

            sesub0 = at(nl)*c10.*(cnew.*clide1);
            sigdif2 = (sigdif1-sesub0);
            DJ(nl+4*NLS) = (sigdif2*sigdif2'-J0)/del;
        end
        figure(3);subplot(2,1,1);title('gradients computed');drawnow;
        % define direction of 1-D search for SD or CG (CG = 1):
        if NNN1 ~= NNN0 
            g=-DJ;gold=g;h=g;
            figure(3);subplot(2,1,1);plot(h);drawnow;
        else
            gold=g;g=-DJ;
            ga=((g-gold)*g')/(gold*gold');
            h=g+ga*h;
            figure(3);subplot(2,1,2);plot(h);drawnow;
        end
        ni=1;nj=0;
        while ni < 5 && increm > 1e-20
            %p1=p;
            ni=1;J1=0;delp=0;Jold=J0;delJ=-1;
            while delJ<0
                delp=delp+ni*increm;de(ni)=delp;
                p1(1:2*NLS)=p(1:2*NLS)+delp*(h(1:2*NLS)+1i.*h(2*NLS+1:4*NLS));
                p3(1:NLS) = p2(1:NLS) + delp*(h((4*NLS)+1:5*NLS));
                sesto=xyvalcs(p1,p3,psis,NNN1,NLS,wfr,wfi,wfr1,time1,time2);  % <<<<< estimate the signal
                sigdif = (FID3(1:NNN1)-sesto);
                J1=sigdif*sigdif';JR(ni)=J1;
                figure(1);subplot(2,1,1);plot((1:ni),JR(1:ni));title(num2str(NNN1));drawnow;
                delJ=J1-Jold;Jold=J1;
                ni=ni+1;nj=nj+1;
            end
            increm=increm/3;
        end
        figure(3);subplot(2,1,1);title('minimum bracketed');pause(0.3);
        R=.61803399;C=1-R;
        if ni >4; 
            x0=de(ni-3);
        else
            x0=0;
        end
        x3=de(ni-1);bx=delp/3;cx=delp;
        x1=bx;x2=bx+C*(cx-bx);
        p1(1:2*NLS)=p(1:2*NLS)+x1*(h(1:2*NLS)+1i.*h(2*NLS+1:4*NLS));
        p3(1:NLS) = p2(1:NLS) + x1*(h((4*NLS)+1:5*NLS));
        sigdif = (FID3(1:NNN1)-xyvalcs(p1,p3,psis,NNN1,NLS,wfr,wfi,wfr1,time1,time2));
        JT1=sigdif*sigdif';
        p1(1:2*NLS)=p(1:2*NLS)+x2*(h(1:2*NLS)+1i.*h(2*NLS+1:4*NLS)); 
        p3(1:NLS) = p2(+1:NLS) + x2*(h((4*NLS)+1:5*NLS));
        sigdif = (FID3(1:NNN1)-xyvalcs(p1,p3,psis,NNN1,NLS,wfr,wfi,wfr1,time1,time2));
        JT2=sigdif*sigdif';
        while abs(x3-x0)>.02*(abs(x3)+abs(x0));
            if JT2<JT1
                x0=x1;x1=x2;x2=R*x1+C*x3;
                JT1=JT2;
                p1(1:2*NLS)=p(1:2*NLS)+x2*(h(1:2*NLS)+1i.*h(2*NLS+1:4*NLS));
                p3(1:NLS) = p2(1:NLS) + x2*(h((4*NLS)+1:5*NLS));
                sigdif = (FID3(1:NNN1)-xyvalcs(p1,p3,psis,NNN1,NLS,wfr,wfi,wfr1,time1,time2));
                JT2=sigdif*sigdif';
            else
                x3=x2;x2=x1;x1=R*x2+C*x0;
                JT2=JT1;
                p1(1:2*NLS)=p(1:2*NLS)+x1*(h(1:2*NLS)+1i.*h(2*NLS+1:4*NLS));
                p3(1:NLS) = p2(1:NLS) + x1*(h((4*NLS)+1:5*NLS));
                sigdif = (FID3(1:NNN1)-xyvalcs(p1,p3,psis,NNN1,NLS,wfr,wfi,wfr1,time1,time2));
                JT1=sigdif*sigdif';
            end
        end
        if JT1 < JT2
            delpf=x1;CC(niter)=JT1;
        else
            delpf=x2;CC(niter)=JT2;
        end
        p(1:2*NLS)=p(1:2*NLS)+delpf*(h(1:2*NLS)+1i.*h(2*NLS+1:4*NLS));
        p2(1:NLS) = p2(1:NLS)+delpf*(h((4*NLS)+1:5*NLS));
        for nl=1:NLS
            if abs(wfr*real(p(nl+NLS))/delt) > 150;
                p(nl)=.05*p(nl);p(nl+NLS)=0;
            end
            if real(p(nl+NLS))>0; 
                p(nl+NLS)=0+1i*imag(p(nl+NLS));
            end
        end
        for nl = 1:NLS
            if abs(wfr*p2(nl)/delt)>150;
                p2(nl)= 0.05*p2(nl);
            end
            if p2(nl)>0
                p2(nl) = 0;
            end
        end
        figure(1);
        subplot(2,1,2);
        plot(log10(CC(1:niter)));title(num2str(CC(niter)));drawnow;
        for nl=1:NLS
            amp(ia(nl),ja(nl))=p(nl);
            expt(ia(nl),ja(nl))=p2(nl);
            expol(ia(nl),ja(nl))=p(nl+NLS);
        end
        ampf(:,:,niter)=fftshift(ifft2(fftshift(WRA.*fftshift(fft2(fftshift(amp))))));
        expf(:,:,niter)=fftshift(ifft2(fftshift(WRA.*fftshift(fft2(fftshift(expol))))));
        exptf(:,:,niter)=fftshift(ifft2(fftshift(WRA.*fftshift(fft2(fftshift(expt))))));
        %     expf(:,:,niter)=expo1;
        %         exptf(:,:,niter)= expt;
        for nl=1:NLS;
            p(nl)=ampf(ia(nl),ja(nl),niter);
            p(nl+NLS) = expf(ia(nl),ja(nl),niter);
            p2(nl) = exptf(ia(nl),ja(nl),niter);
        end      
        %R2f(:,:,niter)=real(-wfr*squeeze(expf(:,:,niter)/delt));
        fr=wfi*imag(squeeze(expf(:,:,niter)))/(twopi*delt);
        figure(4);
        subplot(2,2,1);imagesc(fr);axis image;
        subplot(2,2,2);plot((1:niter),squeeze(imag(wfi*expf(33,34,1:niter)))/(twopi*delt),'k');
        plot((1:niter),squeeze(abs(ampf(27,28,1:niter))),'r--');
        plot((1:niter),squeeze(imag(wfi*expf(27,28,1:niter)))/(twopi*delt),'k');
        hold off;
        subplot(2,2,3);plot(-wfr*real(squeeze(expf(:,33,niter))/delt));title('est. R2*');
        subplot(2,2,4);imagesc(abs(squeeze(ampf(:,:,niter))));axis image;colormap gray
        subplot(2,2,4);title(strcat('estimated magnitude, ',num2str(niter),' iterations'));drawnow;
        figure(5);
        subplot(2,2,1);plot(-wfr1*real(squeeze(exptf(:,33,niter))/delt));title('est. R2');

        NNN0=NNN1;
        CCold=CCnew;
        CCnew=CC(niter);
        niter=niter+1;
    end
end     

R2e=real(-wfr*squeeze(expf(:,:,niter-1)))/delt;% R2' map 
R2p=real(-wfr1*squeeze(exptf(:,:,niter-1)))/delt; %R2 map

frmap=wfi*imag(squeeze(expf(:,:,niter-1)))/(2*pi*delt); %Frequency map
M0=abs(squeeze(ampf(:,:,niter-1)));M0MSK=(M0>.3*max(max(M0)));%Magnitude map
figure(9);subplot(1,2,1);imagesc(M0);axis image
subplot(1,2,2);imagesc(frmap);axis image;colormap gray
subplot(1,3,3);imagesc(R2e);axis image;colormap gray;
figure(10);
subplot(1,2,1);imagesc(R2e);axis image;colormap gray;
subplot(1,2,2);imagesc(R2p);axis image;colormap gray;
toc
   